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introduction 

During  the  interior  ballistic  cycle  in  conventional  tube  weapons,  processes  of 
different  importance  occur.  The  most  important  process,  the  motion  of  the 
projectile,  depends  on  the  burning  of  the  propellant  and  thus  the  generation  of 
gas  pressure  in  the  gun.  The  interior  ballistician  however  knows,  that  beside  this 
key  process  a  series  of  further  processes  must  be  considered,  limiting  the 
design  of  a  charge.  Important  points  of  view  are  for  instance  the  life-time  of  gun 
tubes,  which  is  coupled  to  the  erosion,  or  the  cook-off  of  charges.  These 
examples  show  that  the  generation  of  heat  and  heat  transport  processes  play  an 
important  role.  The  heat  transport  processes  of  heat  conduction,  convection  and 
radiation  as  well  as  heat  transfer  should  be  considered. 

It  is  the  subject  of  this  contract,  to  explore  the  heat  transport  processes  in 
interior  ballistic  flows  by  means  of  numerical  simulation.  The  numerical  simula¬ 
tion  for  high  velocity  flows  under  high  pressure  and  high  temperature  should  be 
limited  however  to  pure  gas  flows  behind  a  projectile.  The  general  subject  of  the 
contract  "Viscous  Modeling  of  the  Interior  Ballistic  Cycle"  emphasizes  the 
viscosity  of  the  interior  ballistics  gas  flow  to  be  included.  The  viscosity  of  the 
flow  is  taken  into  consideration  by  using  the  full  Navier-Stokes  equations  for 
describing  a  multidimensional  (axisymmetric)  unsteady  and  compressible  flow. 
This  kind  of  mathematical  description  opens  the  possibility  of  calculating  the 
formation  of  the  boundary  layer,  which  is  important  for  the  heat  transfer  and  for 
including  turbulence  as  an  essential  process  for  the  heat  convection.  Actually, 
this  kind  of  simulation  considers  the  complete  coupling  of  the  flow  between  the 
core  and  the  transient  boundary  layer  as  well  as  the  interaction  of  the  heat 
transfer  into  the  tube  wall  and  the  wall.  This  report  gives  account  of  the 
achieved  results  of  the  numerical  simulation.  The  comparisons  between  laminar 
and  turbulent  flow  have  been  published  already  at  other  places,  among  others  in 
six  periodic  reports,  and  therefore  are  presented  here  only  in  summary.  The 
investigation  of  the  heat  transfer  to  the  tube  wall,  however,  is  treated  more 
comprehensively. 

Substantial  work  on  this  contract  has  been  performed  by  the  Research  Asso¬ 
ciates 

Dr.  J.  Garloff  (Apr  87  -  March  89)  and 
Dr.  S.  Kocaaydin  (June  90  -  May  91). 
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Numerical  Techniques 

For  the  numerical  solution  of  the  Navier-Stokes  equations  the  DELTA-Code  was 
available.  Originally,  this  computer  code  was  supposed  to  provide  a  numerical 
approximation  of  the  solution  of  the  averaged  two-phase  (gas-solid),  multi¬ 
dimensional  equations  governing  viscous  interior  ballistic  flows  within  conven¬ 
tional  guns  [1,2].  The  numerical  method  is  based  on  a  linearized  Alternating 
Direction  Implicit  (ADI)  scheme  of  the  Briley-McDonald  type  [3]  and  applied  to  a 
moving  coordinate  system  caused  by  the  moving  projectile.  The  implicit  finite 
difference  method  is  very  well  suited,  since  a  satisfactory  resolution  of  the 
boundary  layer  makes  necessary  very  small  grids  in  the  wall  proximity.  To  keep 
the  amount  of  grid  points  as  small  as  possible  a  distribution  of  nonequidistant 
grids  is  used  in  radial  direction  with  most  narrow  grids  at  the  wall  and  increas¬ 
ing  grid  size  in  direction  to  the  tube  axis.  In  axial  direction  equidistant  grids 
were  sufficient.  The  code  was  developed  at  the  BRL  on  a  CDC  mainframe 
computer  in  UPDATE  format.  At  the  time  DELTA  was  transfered  to  the  VAX 
computers  accessible  to  the  Ernst-Mach-lnstitut  (EMI),  only  those  code  modules 
for  an  axisymmetric  laminar  gas  flow  for  adiabatic  walls  had  been  tested.  Beside 
the  main  program  the  graphics  package  based  on  DISSPLA  also  was  installed  on 
the  EMI-computer  for  processing  the  numerical  results.  The  graphic  was 
extended  substantially  to  make  possible  a  sufficient  monitoring  of  the  computer 
runs  and  a  postprocessing. 


Turbulence 

Computational  results  for  the  laminar  interior  ballistic  gas  flow  showed  very  thin 
boundary  layers  [4]  for  the  axial  flow  velocity  as  well  as  for  the  gas  temper¬ 
ature.  By  comparing  the  experimental  results  from  ballistic  simulators  [5,6,7]  to 
the  computational  results  of  the  laminar  flow,  differences  appeared  for  the 
shape  of  the  velocity  boundary  layer  even  for  low  velocities.  The  measured 
velocity  traces  pointed  to  turbulence  as  the  source  of  the  differences. 
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In  a  first  step  we  studied  the  applicability  of  two  algebraic  eddy  viscosity 
models  of  turbulence  using  a  two-layer  approach  since  their  implementation  did 
not  cause  too  many  modifications  in  the  code.  Both  the  Cebeci-Smith-Mosinskis 
model  and  the  Baldwin-Lomax  model  did  not  generate  realistic  turbulent  flow 
profiles  in  the  boundary  layer.  However  the  boundary  layer  thickness  (99  % 
criterion)  looked  similar  to  the  experimental  results.  Although  performing  an 
extensive  parameter  analysis  we  did  not  get  satisfying  results  when  comparing 
with  the  experimental  results  with  the  gun  simulators  of  the  Imperial  College, 
London  [5,6].  In  order  to  describe  adequately  the  heat  transfer  it  is  necessary  to 
start  from  realistic  boundary  layer  profiles.  The  results  using  the  algebraic 
turbulence  models  were  published  in  detail  [8-12]. 

The  experience  with  the  algebraic  models  of  turbulence  caused  us  to  incorpo¬ 
rate  a  more  advanced  turbulence  model  based  on  transport  equations  for  turbu¬ 
lent  quantities.  As  next  step  we  modified  the  widely  used  k-e  model  taking  into 
consideration  compressibility  effects  (131.  Special  emphasis  was  laid  on  the 
modeling  of  pressure-velocity  correlations  appearing  in  the  equation  for  internal 
energy  as  well  as  turbulent  kinetic  energy  k.  This  model  was  strictly  extended  In 
order  to  include  real  gas  effects  (Noble-Abel  equation  of  state)  and  to  facilitate 
the  investigations  on  heat  transfer.  Then  the  k-e  model  was  implemented  in  the 
DELTA  code.  Test  runs  showed  that  the  values  of  the  turbulent  kinetic  dissi¬ 
pation  rate  rose  too  fast  in  the  corner  between  the  projectile  base  and  the  tube 
wall. 

In  view  of  the  numerical  difficulties  with  tho  k-e  model  a  one-equation  model 
was  considered  as  a  more  appropriate  cas .didate.  This  equation  is  a  transport 
equation  for  the  turbulent  kinetic  energy,  .is  dissipation  rate  is  given  by  means 
of  a  mixing  length  specified  algebraically.  By  applying  this  model  in  DELTA  we 
achieved  reasonable  computational  results  in  good  agreement  with  the  experi¬ 
mental  ones.  The  basic  equations  as  well  as  the  comparison  with  the  experi¬ 
ments  are  documented  in  (11,12). 
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Keat  Transfer 

Besides  the  modeling  of  turbulent  gas  flow  close  to  the  wall  the  calculation  of 
the  he£t  transfer  from  the  hdt  propellant  gases  to  the  tube  wall  was  a  further 
important  goal  of  this  contract.  Two  ways  of  solution  were  traced.  First,  the 
boundary  condition  at  the  wail  surface  was  formulated  as  a  heat  balance. 
Second,  a  finite  volume  formulation  was  chosen  such  that  a  volume  across  the 
boundary  Includes  a  part  of  the  flow  and  of  the  wall.  This  method  is  described 
in  detail  by  C.  S.  Kocaaydin  in  the  Appendix  of  this  report  whereas  the  first 
method  is  explained  in  the  following.  Both  methods  required  substantial  exten¬ 
sions  of  the  DELTA  code  ar  well  as  of  the  graphics  for  postprocessing  and 
consequently  plenty  of  coding. 

In  order  to  calculate  the  heat  transfer  to  the  tube  the  boundary  condition  used 
so  far  for  an  adiabatic  wall  surface  was  changed  such  that  the  heat  exchange 
between  flow  and  wall  is  described.  In  this  case  both  the  heat  flux  from  the  hot 
gas  and  the  quantitv  of  radiated  heat  to  the  wall  are  supposed  to  equal  the  heat 
flux  penetrating  into  the  wall: 

3T  3T 

3r  (*^0  ”  3r 

and  on  the  inner  tube  :.urface  the  gas  temperature  Tg  is  supposed  to  equal  the 
wall  temperature  Twi 

Tg  =  T,  . 


The  variables  A,  f,  a  and  r  denote  the  thermal  conductivity,  the  emissivity  of  the 
wall  surface,  the  Stefan-Boltzmann  constant  and  the  radial  coordinate,  respec¬ 
tively.  The  index  g  stands  for  gas,  and  w  for  wall.  The  temperature  Tgo  means 
the  maximum  gas  temperature  in  a  particular  cross-section  of  the  tube.  In  the 
wall  the  two-dimensional  heat  conduction  equation 


K  IF 


is  solved  using  the  same  numerical  scheme  as  for  the  gas  flow  region,  p,  c,  t,  z 
denote  the  density  and  specific  heat  of  the  wall  material,  time  and  axial  coordi¬ 
nate  In  the  wall,  respectively.  At  the  outer  tube  wall  surface  the  simple  engi¬ 
neering  condition 
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is  chosen  where  we  are  using  a  heat  transfer  number  h  and  an  outer  ambient 
temperature  Tamb*  ^  sophisticated  condition,  at  least  for  a  single  shot 
weapon,  is  not  necessary  because  the  heat  usually  does  not  reach  the  outer 
surface  during  the  ballistic  cycle  because  of  the  wall  thickness.  Two  additional 
boundary  conditions  are  needed  in  axial  direction  for  the  heat  conduction  equ¬ 
ation.  At  the  projectile  ba.se,  we  set 


'  \ 


across  the  wall  thickness,  that  is,  we  assume  the  projectile  moves  into  a  cold 
area  with  ambient  temperature.  This  assumption  is  not  true  if  there  is  a  gas 
leakage  around  the  projectile.  At  the  breech,  an  adiabatic  condition 
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seems  to  be  adequate. 


Numerical  Results 
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In  the  following  numerical  results  will  be  presented.  In  order  to  avoid  an  initial 
discontinuity  in  temperatu-'e  at  the  wall  surface  the  initial  temperature  is 
assumed  to  be  ambient  eveiy where  in  the  flow  area  and  in  the  wail.  Also,  the 
initial  pressure  is  ambient.  For  producing  a  time-dependent  pressure  and  temper¬ 
ature  profile  as  in  a  real  tube  weapon,  mass  and  energy  are  added  to  the  one- 
phase  flow  via  source  terms.  These  source  terms  take  into  account  an  e.mptrical 
burning  law  as  a  function  of  pressure.  The  sources  move  with  the  flow.  A  con¬ 
stant  time-step  of  10  //s  was  chosen.  The  computational  mesh  consists  of  97 
uniformly  spaced  grid  points  in  the  axial  direction  and  73  nonuniformly  spaced 
grid  points  in  the  radial  direction  for  the  gas  flow  whereas  73  uniformly  spaced 
grid  points  are  used  in  the  wall  In  radial  direction. 


Table  1 .  Parameters  of  the  Test  Run 


Bore  Diameter 

20  mm 

Tube  Length  (Breech  to  Muzzle) 

2.0  m 

Initial  Chamber  Length 

0.175  m 

Projectile  Mass 

120  g 

Ratio  of  Specific  Heats 

1.271 

Covoiume 

1.08  •  10-3  m3/kg 

Molar  Mass 

23.8  g/mole 

Barrel  Material 

Steel  4130 

Density 

7800  kg/m3 

Thermal  Conductivity 

43  W/m/K 

Specific  heat 

460  J/kg/K 

Emissivity  Factor 

0.7 

Initial  Gas  Pressure 

0.1  MPa 

Initial  Gas  and  Wall  Temperature 

293  K 

Initial  Gas  Velocity 

0  m/s 

Shot  Start  Pressure 

0.1  MPa 

The  smallest  grid  size  in  radial  direction  is  0.61  //m  on  the  gas  side  and  34.7  ;/m 
in  the  wall.  Other  parameters  for  this  test  run  are  summarized  in  Table  1. 
Several  Figures  referring  to  the  boundary  layer,  wall  temperature  distribution, 
heat  transfer,  spatial  distribution  of  flow  quantities  and  projectile  motion  will  be 
presented.  5.81  ms  after  the  beginning  of  the  shot  cycle  the  projectile  leaves 
the  muzzle  with  the  velocity  of  941  m/s  having  run  a  distance  of  1 .825  m  down 
the  tube.  Bore  friction  is  neglected. 

Figures  1  to  3  show  information  on  the  projectile  motion.  Pressure  as  well  as 
temperature  rise  and  drop  versus  time  are  displayed  in  Figures  4  and  5,  respec¬ 
tively.  As  mentioned  before  mass  and  energy  is  added  to  the  flow  in  order  to 
achieve  realistic  gas  and  temperature  profiles.  This  is  based  on  a  hypothetical 
mass  of  propellant  the  consumption  of  which  is  shown  in  Figure  6.  The  spatial 
distribution  of  the  axial  flow  velocity,  radial  flow  velocity  and  gas  pressure  are 
plotted  in  Figures  7  to  9,  respectively,  for  the  time  of  muzzle  clearance.  At  this 
time  the  boundary  layer  in  the  front  part  behind  the  projectile  is  fully  developed 
(cf.  Fig.  7).  The  turbulent  kinetic  energy  and  the  Reynolds  stress  are  shown  in 
Figures  10  and  11,  respectively. 
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In  Figure  1 2  the  temperature  distribution  in  the  gas  and  in  the  wall  are  given  for 
the  same  time.  The  wail  surface  is  the  line  where  the  red  and  green  areas  meet 
each  other.  The  wall  temperature  is  resolved  more  detailed  in  Figure  13. 
Temperature  measurements  are  taken  at  fixed  locations  along  the  wall.  There¬ 
fore,  the  Figures  14  to  16  show  at  the  positions  100,  200  and  300  mm  off  the 
breech  the  temperature  curves  for  the  wail  surface  and  for  the  depths  in  the 
wall  of  10,  30,  so,  100 //m.  One  has  to  keep  in  mind  that  the  initial  position  of 
the  projectile  is  175  mm  off  the  breech.  The  wall  temperature  depends  on  the 
heat  flux.  The  heat,  transfered  from  the  hot  gas  to  the  cool  wall  along  the  tube 
is  integrated  over  time  in  Figure  17.  The  decrease  of  the  transfered  heat  from 
the  maximum  at  about  200  mm  is  caused  by  the  decrease  in  heat  convection 
when  approaching  to  the  breech.  The  convection  of  heat  is  supposed  to  be  the 
dominant  heat  transport  process  to  the  wall  in  order  to  feed  the  gas  layer  next 
to  the  wall  with  heat.  This  is  one  of  the  most  interesting  results  of  this  heat 
transfer  study. 

Several  additional  figures  illustrate  some  other  features  of  the  viscous  flow 
solution.  Next,  the  development  of  the  velocity  boundary  layer  (99  %  criterion) 
is  shown  in  a  sequence  of  plots  in  Figures  18  to  22.  For  the  time  of  muzzle 
clearance  the  radial  profiles  of  the  axial  flow  velocity,  radial  flow  velocity  and 
gas  and  wall  temperature  are  displayed  occuring  0.5,  1.5  and  1.95  m  off  the 
breech.  This  concludes  the  presentation  of  numerical  results  fc  the  particular 
20-mm  test  example.  Of  course  there  exist  many  other  interesting  results.  The 
figures  and  explanations  given  here  are  supposed  to  represent  the  most 
important  attributes  of  the  viscous  turbulent  interior  ballistics  one-phase  flow 
with  heat  transfer. 

As  another  test  example  the  22.2  mm  experimental  gun  simulator  of  the 
Imperial  College  was  chosen.  The  essential  parameters  of  this  test  rig  are  given 
in  Table  2.  The  gas  used  is  low  grade  nitrogen.  The  flow  is  a  pure  expansion 
flow  starting  from  the  initial  gas  pressure  of  1.1  MPa  behind  the  projectile.  This 
example  has  been  used  for  an  extensive  comparison  of  the  particular  flow 
characteristics  between  experiment  and  numerical  simulation.  Above  all  the 
modeling  of  the  turbulence  was  verified.  Since  those  results  have  been 
published  in  Refs.  10  to  12  only  the  heat  transfer  results  are  presented  in  the 
following.  They  are  refered  to  the  Figures  32  to  41.  The  computational  mesh 
again  uses  97  x  73  grid  points.  Figures  32  to  34  show  the  projectile  motion. 
Because  of  the  low  initial  pressure  the  achieved  maximum  velocity  is  only  about 
200  m/s. 
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Table  2.  Parameters  of  the  gun  simulator 


Bore  Diameter  22.2  mm 

Tube  Length  (Breech  to  Muzzle)  3.56  m 

Initial  Chamber  Length  1 .3  m 

Projectile  Mass  10.5  g 

Ratio  of  Specific  Heats  1 .41 6 

Covolume  0  m^/kg 

Molar  Mass  28.0  g/mole 


Barrel  Material  and  Properties  as  in  Table  1 


Initial  Gas  Pressure  1.1  MPa 

Initial  Gas  and  Wall  Temperature  293  K 

Initial  Gas  Velocity  0  m/s 

Shot  Start  Pressure  0.1  MPa 


The  pressure  as  well  as  the  temperature  drop  both  at  the  projectile  base  and 
breech  according  to  the  expansion  are  plotted  in  the  Figures  35  and  36.  Since 
the  expansion  starts  from  ambient  temperature  the  gas  temperature  drops  below 
the  ambient  temperature.  This  means  heat  is  transfered  from  the  wall  to  the 
flow  as  shown  in  the  Figures  37  to  40  where  the  wall  temperature  at  the  loca¬ 
tions  1.0,  1.2,  1,5  and  2.0  off  the  breech  are  shown  for  the  wall  surface  and 
the  wall  depths  of  10,  30,  50  and  100  //m.  The  integral  heat  flux  (cf.  Fig.  41), 
therefore,  is  negative.  The  computed  temperature  drops  in  the  wall,  however,  is 
not  in  agreement  with  those  experimental  results  published  in  Ref.  17].  There, 
the  measured  temperature  decreases  about  4  to  6  °C  whereas  the  computed 
temperature  decreases  about  1  °C  or  less.  The  reasons  for  the  difference  are 
not  known  at  this  writing. 
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Conclusions 

The  Navier-Stokes  solution  is  able  to  predict  particular  flow  features  of  a  one- 
phase  interior  ballistic  flow  such  as  boundary  layer,  turbulence  or  heat  trans¬ 
port.  A  one-equation  model  of  turbulence  seems  to  be  best  suited  to  incorporate 
turbulence  important  to  generate  the  proper  amount  of  heat  convection. 
Another  satisfying  result  concerns  the  success  in  calculating  the  rise  of  temper¬ 
ature  in  the  tube  wall. 

Some  of  the  numerical  results  need  to  be  verified  experimentally  above  all  in  the 
high  velocity  regime.  An  extension  of  the  study  to  the  application  of  this  mathe¬ 
matical  model  on  two-phase  flows  with  inert  and  reactive  solid  particles  is 
recommended. 
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Axial  flow  velocity  at  muzzle  clearance 
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Fig.  8:  Radial  flow  velocity  at  muzzle  clearance 


TIME-STEP=  582  TIME=  5.809[ms] 
PROJ.DISPLACEMENT=2.000[m]  PROJ.VEL.=  941[iVs] 
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Rg.  9: 


Gas  pressure  at  muzzle  clearance 


TIME-STEP=  582  TIME=  5.809[ms] 

PROJ .DISPLACEMENT=2.000[m]  PROJ.VEL.=  94l[n/s] 
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Fig.  10:  Turbuler<t  kinetic  energy  at  muzzle  clearance 


TIME-STEP=  582  TIME=  5.809{ins] 
PROJ.DISPLACEMENT=2.000[in]  PROJ.VEL.=  941  [11/3] 
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Fig.  11: 


Reynolds  stress  at  muzzle  clearance 
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Spatial  tamparatura  distribution  in  gat  and  waii 
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Fig.  13: 


Wall  temperature  at  muzzle  clearance 
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Fig.  14:  Temperature  In  the  wall  versus  time  100  mm  off  the  breech  for  the 
surface  and  the  depths  10,  30,  50  and  lOOivm 
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Fig.  15: 


Temperature  In  the  wall  versus  time  200  mm  off 
surface  and  the  depths  10,  30,  50  and  100//m 
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Fig.  18:  Velocity  boundary  layer  at  3.9  ms 
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Fig.  19:  Velocity  boundary  layer  at  4.4  ms 


TIME-STEP=  490 
PROJ.  POSITION=1.187[m], 


[raui] 

ssaNxoiHL  H3AV1  AHvamoa 


Fig.  20:  Velocity  boundary  layer  at  4.9  ms 
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g.  21:  Velocity  boundary  layer  at  5.4  ms 
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Fig.  23:  Radial  profile  of  the  axial  flow  velocity  0.5  m  off  the  breech  at 
muzzle  clearance 
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Fig.  24:  Radial  profile  of  the  axial  flow  velocity  1.5  m  off  the  breech  at 
muzzle  clearance 
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Fig.  25:  Radial  profile  of  the  axial  flow  velocity  1.95  m  off  the  breech  at 
muzzle  clearance 
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Fig.  27:  Radial  profile  of  the  radial  flow  velocity  1.5  m  off  the  breech  at 
muzzle  clearance 
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Fig.  28;  Radial  profile  of  the  radial  flow  velocity  1 .95  m  off  the  breech  at 
muzzle  clearance 
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Radial  profile  of  the  temperature  in  gas  and  wall  0.5  m  off  the 
breech  at  muzzle  clearance 
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Radial  profile  of  the  temperature  In  gas  and  wall  1  95  m  nf#  tho 
breech  at  muzzle  clearance 


TIME  HISTORY 


PROJ.  TRAVEL  [m] 


TIME  HISTORY 


TIME  [m; 


TIME  HISTORY 


TIME  HISTORY 


-  53- 


TIME  [m; 


[ME 


L6.360[ms] 
EL=  202.3rm 


-  57 


APPENDIX 


INTERIOR  BALLISTICS 
AND 

GUN  BARREL  HEAT  TRANSFER 


Development  of  DELTA-Code  considering 
Unsteady  Heat  Conduction 


Dr.-Ing  C.  S.  Kocaaydin 


9.10.1991 


-  58  - 


1.  Outline 

The  purpose  of  this  research  project  was  to  develop  a  gun  barrel  heat  transfer 
model  in  order  to  predict  the  thermal  boundary  layer  development  and  the 
heat  load  to  the  barrel.  This  problem  was  actualley  embedded  into  the  deve¬ 
lopment  of  a  general  computer  code  in  order  to  analyse  the  overall  interior 
ballistics  cycle  in  detail. 

During  the  period  from  June  90  to  November  90,  the  major  research  activity 
was  concentrated  on  the  development  of  a  new  method  to  consider  the 
heat-transfer  to  the  gun  wall.  The  theoretical  aspects  of  this  new  method 
have  been  documented  already  in  the  last  report  /!/.  In  this  final  report,  we 
are  going  to  summarize  the  progress  achieved  in  the  research  project  during 
the  period  from  December  90  to  September  91. 

2.  Experience  with  the  DELTA  program 

The  development  of  a  gun-barrel  heat  transfer  model  and  its  implementation 
into  the  DELTA  code  encountered  several  difficulties  due  to  the  numerical 
stability  problems.  This  situation  has  led  to  some  work,  which  was  initially 
not  considered. 

2.1  Numerical  Instability 

Due  to  the  central  differences  and  excessive  quasi-linearization  used  within 
the  discretization  methodology  of  the  DELTA  code,  this  program  cannot  be 
considered  as  robust.  Therefore,  it  is  always  necessary  to  adapt  carefully 
the  grid  and  the  time-step  for  every  initial  condition.  Otherwise,  the  results 
are  suffering  from  some  high-frequency  oscillations.  On  the  hand,  to  combine 
the  DELTA  program  with  a  new  boundary  algorithm  has  made  the  stability 
range  narrower.  The  computation  of  thermal  boundary  layer  leads  to  high 
gradients  in  temperature  distribution  within  the  computation  domain.  So  we 
have  encountered  some  basic  problems  during  this  research  period. 

As  it  is  well  known,  the  remedy  to  the  instability  is  some  sort  of  artificial 
viscosity.  Due  to  the  basic  formulation,  which  is  currently  applied  in  the 
DELTA  code,  it  is  necessary  to  put  explicit  and  implicit  artificial  viscosity. 
Although  the  range  of  explicit  viscosity  terms  are  bounded  upwards,  there 
is  actually  no  limit  given  to  the  implicit  artificial  terms.  After  several  tests 
have  been  run  with  the  DELTA  code,  it  was  confirmed,  that  the  DELTA  code 
doesn't  react  stable  if  the  coefficient  of  implicit  artificial  viscosity  is 
increased.  It  is  recommended  to  consider  some  recent  publications  on  the 
role  of  artificial  viscosity  /2/. 

Af  ter  a  rigorous  control  of  the  program,  it  was  found  that  this  above  mentioned 
phenomenon  is  a  result  of  an  error  in  the  implementation  of  implicit  artificial 
viscosity.  Let's  consider  the  linearized  block  implicit  scheme  used  in  the 
DELTA  code.  Complete  description  of  the  method  can  be  found  in  /3,4/. 
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The  general  form  of  the  governing  nonlinear  partial  differential  equation 
is 


^  *  G(r.  t.  y.  y^. 


y.*  y 


rr 


rz 


(1) 


where  the  independent  variables  of  time,  radial  and  axial  coordinates  are 
denoted  by  t,  r,  z,  respectively.  The  vector  of  dependent  variables  is 
denoted  by  y.  The  alternating  direction  -  implicit  method  (ADI)  used  in 
DELTA  is 


[l  -  3  D”]y*  =  [l  -  3D^’]  y^  ♦  Ly^ 

[l  -  0  D;]y'*=  [l  -  y^  ♦{y*- y*^)  .  (2) 

The  operator  matrixes  D^.  and  have  the  following  structure 
D^=  (D.  DR,  DRR)  , 

(DZ.  DZZ)  .  (3) 


with  the  Jacobian  matrices 


d  =  4^ 

dy 


DR  = 


ciG 


DRR 


dG 

^y. 


rr 


D2  = 

dy 

•^Z 


DZZ  = 


iG 

<^y,. 


In  the  present  code,  the  implicit  artificial  viscosity  terms  were  added  to  the 
operator  matices  DRR  and  DZZ  on  each  side  of  the  equations  -  see  subroutines 
SWRINT  and  SVVZINT  But  it  is  only  allowed  to  use  these  terms  on  the  left 
implicit  side. 


After  this  error  has  been  eliminated,  the  DELTA  code  has  become  really 
stable.  It  should  be  clear,  that  the  use  of  artificial  viscosity  leads  to  non¬ 
physical  diffusion.  Therefore,  these  coefficients  should  be  kept  as  low  as 
possible. 


There  is  another  important  question  within  the  context  of  the  stability 
problem.  This  is,  how  the  Jacobian  matrices  should  be  computed.  It  Is  possible 
to  compute  them  analytically  or  numerically.  A  definite  answer  to  this 
question  cannot  be  given  at  this  moment. 


2.2  Reduction  of  computation  time 


In  order  to  reduce  the  computation  time,  it  was  decided  to  use  only  exact 
Jacobi  matrices.  So  the  algorithm  was  completely  reconsidered.  The  compu¬ 
tation  of  the  right-hand-side  and  of  the  Jacobians  has  been  recoded.  By 
storing  several  intermediate  quantities  the  computation  time  has  been 
reduced  drastically.  On  the  other  hand,  due  to  the  higher  robustness  of  the 
code  after  the  correction  of  the  artificial  viscosity,  the  time-steps  could 
also  be  increased.  Some  existing  subroutines  were  further  developed  to 
compute  the  Jacobians  analytically  even  in  case  of  turbulent  flow.  Unfor- 
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3.  Development  of  e  new  boundary  condition  to  couple  the  fluid  flow  and 
heat  conduction  within  the  wall 


In  this  section,  some  details  of  the  formulation  will  be  presented.  These 
are  the  choice  of  the  grid,  numerical  aspects  and  boundary  conditions 

3.1  Choice  of  the  grid  *■  Stability  of  heat-transfer 

The  most  important  aspect  in  the  implementation  of  a  coupled  heat-transfer 
model  was  the  question  concerning  the  grid  resolution  in  the  gas  and  the 
wall.  There  were  some  open  questions.  The  choice  of  any  grid  can  affect  the 
results.  This  uncertainty  has  been  resolved  in  the  following  way 

Actually,  it  is  well  known,  that  the  grid  should  be  fine  enough  within  the 
boundary  layer.  Here,  we  are  considering  how  the  coupling  of  fluid  and  solid 
can  influence  the  unsteady  heat-transfer. 


As  we  are  primarily  interested  in  the  unsteady  heat  transfer,  we  should 
look  carefully  to  the  unsteady  heat  conduction  equation,  where  temperature, 
density,  specific  heat  capacity  and  thermal  conductivity  are  denoted  by  T, 
p,  c  and  k  respectively. 


dt 


with  thermal  diffusivity  a,  defined  as 
k 

a  = - 

pc 

The  heat  equation  can  be  put  into  a  dimensionless  form  through  the  choice 
of  some  appropriate  reference  values. 

d±  _  (5^9 
dt  ”  di)^  ’ 

with 


There  is  an  order  of  magnitude  difference  in  the  diffusivity  of  steel  and  gas. 
So  we  should  also  analyse  the  unsteady  thermal  contact  problem. 


^1  = 
dt 


d^T;, 
2  ax2 


“l  * 


If  the  heat  equation  should  be  resolved  in  two  different  media  simultaneously, 
the  characteristic  length  of  the  domain  should  be  chosen,  such  that  the 
time-scale  of  the  dimensionless  equations  should  be  the  same. 
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This  leads  to  the  following  criteria 

>/?  . 


So  the  computational  domain  used  to  determine  the  boundary  temperature 
must  be  defined  in  a  similar  manner.  Actually,  this  fact  is  also  prescribed  by 
the  stability  criteria  of  an  explicit  integration  scheme.  This  aspect  will  be 
mentioned  in  the  next  section. 

3.2  Numerical  aspects  of  the  formulation 

The  underlying  theory  for  the  formulation  of  a  special  boundary  condition 
to  couple  the  flow-field  inside  the  gun  and  the  heat-transfer  in  the  gun- 
barrel  has  been  previously  reported  /!/.  But  this  numerical  implementation 
has  encountered  several  difficulties.  In  order  to  remove  this  problems,  all 
numerical  aspects  have  been  reconsidered  during  the  second  half  of  the  last 
period.  Here  we  are  giving  a  full  description  of  the  method. 

3.2.1.  Discretisation  and  integration  of  the  energy  equation 

To  compute  the  temperature  at  the  gun-wall,  the  energy  equation  should 
be  integrated  in  a  time-accurate  manner.  The  computation  domain  covers 
naturally  the  gas  flow  inside  the  tube  and  the  heat  conduction  in  the 
gun-barrel.  In  order  to  localize  the  problem,  it  will  be  appropriate  to 
make  some  assumptions.  As  long  as  the  radial  heat-transfer  is  stronger 
than  the  axial  processes,  we  can  ignore  the  axial  convective  and  heat-flux 
terms  in  the  energy  equation.  This  assumption  is  true  in  case  of  high 
Peclet  numbers..  Reducing  the  two-dimensional  problem  into  a  radial  heat- 
transfer  problem,  we  can  sweep  along  radial  direction  independently  at 
each  axial  section. 


Let’s  state  the  problem  once  again.  We  are  looking  to  a  radial  heat-transfer 
problem  and  seeking  only  the  temperature  at  a  single  point,  which  is  the 
common  bounda'y  of  two  different  media.  Choosing  an  explicit  scheme  to 
integrate  the  energy  equation,  the  computational  domain  will  cover  only  a 
small  region  around  the  boundary  point. 

It  should  be  remembered  that  an  explicit  scheme  should  obey  some  stability 
criteria.  In  this  case  the  time  step  has  been  already  prescribed  through  the 
solver  for  the  fluid- dynamics  problem.  So  we  have  to  define  the  integration 
domain  around  the  boundary  point  carefully.  It  is  well  known,  that  the 
stability  criteria  of  an  explicit  scheme  for  a  diffusion  equation  is  given  like 

• 

4x2  2 
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So  we  have  to  consider  the  gas  and  wall  sides  independently.  Due  to  the 
different  material  constants,  there  is  an  order  of  magnitude  difference  in 
the  heat  diffusion  coefficient.  This  is  an  essential  point  for  the  coupling  of 
two  media.  It  should  be  remembered  that  the  diffusion  coefficient  depends 
also  on  temperature.  That's  why  it  is  appropriate  to  use  an  auxiliary  grid 
within  this  boundary  algorithm.  Otherwise,  it  is  a  hopeless  task  to  find  an 
appropriate  grid  to  fulfill  all  of  the  requirements  before  running  the  program. 

Actually  we  are  going  to  formulate  a  three-point  explicit  scheme  for  the 
integration  of  the  energy  equation,  which  is  discretized  according  to  finite- 
volume  method.  Considering  the  same  stability-criteria,  we  compute  the 
location  of  a  northern  node  N  in  the  wall  and  of  a  southern  node  S  in  the 
gas.  The  boundary  node  is  described  by  the  letter  P.  The  next  step  is  the 
computation  of  all  state  variables  at  this  new  defined  points.  A  simple 
interpolation  is  good  enough  for  our  purposes  here. 

ArpN  =  'Jl  At 

Afpg  ='/2Ato<g^5 

The  radial  coordinate  Cp  of  the  boundary  point  P  is  defined  by  the  geometry 
of  the  gun  barrel. 

=  Tp  +  Arpfj  . 

Tg  =  Tp  -  Arps  . 

Due  to  the  difference  in  material  constants,  we  cannot  conside.”  a  single 
composite-cell  at  the  boundary.  Therefore  we  should  define  a  boundary  cell 
for  each  medium.  The  upper  cell  U  is  the  boundary  cell  in  the  gun-wall.  The 
lower  cell  L  is  just  the  opposite,  boundary  cell  in  the  gas.  An  important 
aspect  of  the  chosen  finite -volume  formulation  is  that  we  are  going  to 
make  an  energy  balance  for  each  boundary  cell.  After  computing  the  tempe¬ 
rature  of  cell  centers  at  the  new  time  level,  the  temperature  of  the  boundary 
point  will  be  determined  assuming  an  equilibrium  of  heat  fluxes. 

The  temperature  at  each  boundary  cell  T^  or  T^  depends  directly  on  the 
boundary  temperature  Tp  ,  but  they  are  not  the  same.  At  every  time-step, 
the  state  variables  at  the  center  of  each  boundary  ceil  is  defined  through  a 
weighted  averaging  of  the  adjacent  nodes.  The  weighting  factors  are  given 
just  through  geometric  considerations.  If  we  define  a  cell  center  for  each 
boundary  cell,  so  the  location  can  be  computed  simply  as  given  below. 


!  ! 
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According  to  this  definition,  it  is  also  a  simple  arithmetic  to  define  the 
weighting  factors. 


Let  us  summarize  the  discretization  of  energy  equation  in  radial  direction 
at  each  axial  section.  We  are  considering  now  only  two  cells  U  and  L  with 
temperatures  and  Tj^.  The  simplified  energy  equation  is  given  below, 
using  the  common  notation  for  finite-volume  methods. 


^  [(PT)”*' 


■In  ■  ‘•s.  ■  [{P'r'')N  ■  (P%'')s.  * 


The  discretized  equations  for  a  radial  heat-transfer  sweep  at  an  axial 
section  are  given  below  using  the  same  terminology.  Temperature  T  and 
density  p  will  be  computed  at  the  new  time  level  t”**  .  So  they  are 
unknowns  in  this  formulation.  All  other  variables  as  v,  p,  Cy  will  be  taken 
from  the  latest  available  time  level  t”. 

Energy  balance  in  the  upper  boundary  cell  U: 

Assuming  the  density  of  the  gun  wall  is  constant,  the  only  unknown  at  the 
upper  cell  is  temperature. 


Energy  balance  in  the  lower  boundary  cell  L: 

In  the  gas,  the  density  depends  on  temperature  and  pressure.  We  assume 
that  the  pressure  changes  can  be  neglected  for  a  single  time-step.  So  the 
density  will  be  a  function  of  temperature.  The  following  state-law  for  a 
Nobel-Abel  type  gas  should  be  considered. 


p  =  irr 

MP 


♦  I) 


^  [(pT)r’  -  (pT);:]-  [  -  (v  J  „Fps(pT)"„  ] 

-t;  ))]-[-WpsFpsPps]*®L  \ 


-  65  - 


( 


This  nonlinear  equation  can  be  easily  solved  by  a  damped  Newton  method. 


Finally  we  should  define  the  heat-flux  at  the  boundary  qp,  so  that  the 
equations  can  be  integrated  explicitly.  For  this  purpose,  we  assume  a  local 
equilibrium  of  heat-flux  between  the  wall  and  gas.  Using  the  same  termino¬ 
logy,  this  means  that  the  heat-flux  from  the  cell  center  L  towards  the 
boundary  point  P  is  equal  to  the  heat-flux  from  the  boundary  point  P 
towards  the  upper  cell  center  U.  It  is  easy  to  define  the  heat-flux  in  the 
solid,  because  heat  conduction  is  the  only  heat-tn  nsfer  mechanism.  On 
the  other  hand,  it  is  very  difficult  to  describe  the  heat-flux  through  a  fluid 
layer  at  the  wall.  Considering  the  following  relation  according  to  the  local 
equilibrium 


the  heat  flux  is  defined  simply  using  the  information  on  the  solid  side 


3.2.2  Determination  of  surface  temperature 

After  the  equations  have  been  solved  for  temperature  of  cell  centers  at  the 
new  time  level,  the  final  step  will  be  to  determine  the  new  temperature  at 
the  boundary.  We  have  already  assumed  that  there  is  an  equilibrium  of  heat- 
fluxes  within  each  medium  at  the  boundary  Since  we  have  considered  the 
unsteadiness  of  the  problem  and  also  taken  account  of  all  the  convective 
terms,  this  assumption  can  be  justified. 


In  order  to  make  the  details  clear,  we  should  look  at  first  to  a  simple  contact 
problem.  Let's  consider  the  contact  of  two  solid  layers  (  upper  and  lower  ) 
having  different  heat-conductivity  kpj,^  and  kpg  respectively.  The  thickness 
of  the  layers  will  be  described  by  Arpj^^  and  Arpg.  In  this  case,  heat  conduc¬ 
tion  is  the  only  heat-transfer  mechanism.  Finally  this  composite  layer  should 
be  held  between  a  temperature  difference  T^  and  T^.  The  boundary  tempera¬ 
ture  between  two  .solid  layers  Ci*.~<  be  computed  considering  an  equilibrium 
of  heat-flux  using  the  following  equation 


p 


jn 

^»‘PN 

*^PN 


jn*l^ 


PS 

^ps 

kps 


Actually  this  relationship  cannot  be  used  in  case  of  contact  between  a  solid 
and  fluid  layer.  The  reason  is  that  the  heat-flux  on  the  fluid  side  is  not  only 
due  to  the  heat-conduction.  However  we  can  use  the  already  available  infor¬ 
mation  on  the  velocity  distribution  near  to  the  wall  in  order  to  define  the 
wall  heat-flux  on  the  fluid  side.  This  is  to  say  that  the  method  of  solution 


is  an  analogy  between  the  velocity  and  temperature  distribution  sometimes 
known  as  the  Reynolds-analogy.  In  case  of  high-speed  flows  having  a  Prandtl 
number  close  to  one,  this  assumed  similarity  between  velocity  and  tempera¬ 
ture  distribution  Is  reasonable. 


At  this  point  it  is  important  to  note  that  we  are  going  to  use  the  Reynolds- 
analogy  only  in  the  boundary  cell  which  we  have  already  discussed.  The 
velocity  gradient  at  the  wall  is  evaluated  on  the  original  grid. 
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So  we  can  define  a  heat-transfer  coefficient  a 
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Finally  the  wall  temperature  can  ’  e  computed  similarly. 


du 

dr 
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3.2.3  Influence  of  grid-motion  on  the  wall  temperature 

Up  to  now,  we  have  considered  the  heat-transfer  within  a  stationary  grid 
system.  Therefore  the  wall  temperature  should  be  modified  concerning  the 
grid  motion.  This  is  accomplished  through  a  simple  interpolation  of  the  wall 
temperature  using  the  old  grid  and  computing  the  temperature  at  the  new 
grid  points  along  the  wall.  Concerning  the  wall/base  corner  point  special 
care  is  necessary.  Due  to  the  projectile  motion,  the  grid  is  expanding.  The 
regions  which  are  entering  the  computational  domain  are  assumed  to  be  at 
the  room  temperature.  But  this  is  not  valid  for  the  wall  points,  except  the 
wall/base  corner  point.  Due  to  the  strong  expansion  of  the  grid,  it  can 
happen  that  some  wall  points  will  be  between  the  old  and  new  projectile 
locations.  Therefore  we  have  to  compute  the  temperature  of  the  wall/base 
corner  point  belonging  to  the  old  grid  at  the  new  time  level.  This  value  is 
used  here  only  for  internal  purposes  during  interpolation  of  wall 
temperature.  It  will  be  helpful  to  compare  with  the  details  of  boundary 
conditions  in  section  3.4  . 
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3.3  Diicuaaion  of  the  new  boundary  algorithm. 

The  presented  formulation,  has  the  aim  to  consider  the  unsteady  and 
convective  nature  of  the  heat-transfer  problem  at  the  boundary  of  a  fluid 
and  solid.  It  has  been  observed,  that  a  further  assumption  in  order  to 
simplify  the  problem  can  be  made.  This  is  that  the  radial  convection  and 
the  pressure  work  terms  can  be  neglected,  if  the  radial  velocity  is  truly 
vanishing.  This  is  not  always  the  case.  Even  some  very  small  radial  velocity 
leads  to  very  high  radial  convection  and  pressure  work,  because  both  of 
the  state  variables  temperature  and  pressure  attain  very  high  magnitudes. 
Within  the  context  of  stability,  this  phenomenon  causes  also  trouble. 
Remember  that  we  are  using  only  a  one-step  explicit  scheme.  So  it  may  be 
better,  to  use  a  fine  grid  at  the  boundary  and  neglect  these  terms  as  long 
as  an  explicit  scheme  is  used. 

So  the  simplified  energy  equation  describes  the  unsteady  temperature 
distributior.  in  a  Coutte-flow,  balancing  the  heat  conduction  with  the 
heat-generation  through  viscous  dissipation.  Actually  all  universal  wall- 
functions  have  originated  from  such  models. 

3.4  Modification  of  Boundary  Conditions 

The  definition  of  boundary  conditions  is  usually  considered  as  the  most 
important  aspect  of  a  numerical  simulation  program.  Therefore,  special 
attention  is  required. 

In  order  to  make  the  DELTA  code  transparent,  the  subroutines  concerning 
boundary  conditions  have  been  simplified.  Concerning  the  momentum 
equation,  no-slip  condition  is  used  everywhere  except  the  base,  where  the 
velocity  is  prescribed.  At  the  breech  the  axial  pressure  gradient  and  along 
the  wall  the  radial  pressure  gradient  is  set  to  zero.  Along  the  axis  the 
symmetry  condition  is  well  defined. 

For  the  gasdynamical  problem  in  an  interior  ballistics  cycle  with  a  coupled 
gun-barrel  heat-transfer  modelling,  the  boundary  conditions  at  the  base  of 
the  projectile  were  modified.  It  has  been  observed  that  the  thickness  of  the 
thermal  boundary-layer  is  unrealistically  big  behind  the  projectile  base. 
Therefore  a  correction  of  pressure  and  entropy  is  necessary.  The  flow  at  the 
base  has  been  assumed  as  one-dimensional.  After  every  time-step,  the 
pressure  and  entropy  along  the  base  of  the  projectile  were  set  equal  to 
corresponding  values'at  the  axis.  It  should  be  mentioned,  that  the  grid  is 
uniform  in  the  axial  direction.  That's  why  the  resolution  of  any  viscous 
flow  phenomena  at  the  base  cannot  be  achieved. 

One  important  detail  belongs  to  the  handling  of  the  wall/base  corner  point. 
Actually,  there  are  at  least  two  different  states  at  this  singular  point.  If  we 
consider,  that  this  point  belongs  to  the  wall,  its  velocity  is  set  to  zero. 
Concerning  the  entropy,  the  situation  is  some  what  different.  This  singular 
point  is  also  the  edge  of  the  boundary-layer.  If  we  as.sume  that  the  boundary 
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layer  has  not  been  developed,  then  we  can  define  the  temperature  of  this  point 
through  extrapolation  from  the  gas  flow.  Os  the  other  hand,  if  we  are  going 
to  compute  the  heat  conduction  within  the  gun-barrel,  we  should  remember 
that  the  heat-flux  has  not  entered  into  the  wall  yet.  So  we  should  consider 
that  this  singular  point  is  still  at  room  temperature. 

4.0  Dlscuaalott  of  numerical  results  and  conclualon 

The  gun-barrel  heat-transfer  is  controlled  through  convective  and  viscous 
processes  within  boundary  layer.  At  the  beginning  of  the  ballistic  cycle,  the 
wall  temperature  close  to  the  base  of  the  projectile  attains  very  high  values. 
This  peak  is  caused  through  the  developing  boundary  layer.  At  the  edge  of 
the  boundary-layer  the  viscous  dissipation  is  very  high.  Beside  that  the  hot 
gas  of  the  core-flow  touches  instantanousiy  the  wall/base  corner  point. 

Another  important  process  is  the  radial  convection  behind  the  projectile 
base.  The  axial  velocity  of  the  gas  behind  the  projectile  is  larger  than  the 
projectile  velocity.  Due  to  the  deceleration  of  the  gas  flow  behind  the 
projectile,  there  is  a  change  of  flow  direction  from  axial  towards  radial.  So 
the  heat  load  is  further  increased  through  radial  convection.  It  should  be 
remembered  this  deceleration  process  is  also  an  adiabatic  heating.  The  flow- 
field  behind  the  projectile  is  a  combination  of  two  completely  different 
phenomena.  The  adiabatic  heating  Is  reacting  against  the  expansion  flow 
due  projectile  motion. 

Concerning  the  numerical  simulation  of  interior  ballistics,  there  is  stiil  some 
need  for  the  development  of  an  efficient  method.  The  DELTA  code  is  actually 
describing  the  hydrodynamics  of  the  problem.  Even  for  this  problem,  it  will 
be  appropriate  to  rewrite  the  code.  Using  the  same  theoretical  modelling  with 
a  modern  implicit  discretisation  method,  based  on  upwind  finite  differences 
and  considering  an  efficient  code-structure,  thi*5  task  can  be  fully  achieved 
in  a  reasonable  short  time. 
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